Transport Coefficients for Granular Media from Molecular Dynamics Simulations* 
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Under many conditions, macroscopic grains flow like a fluid; kinetic theory predicts continuum 
equations of motion for this granular fluid. In order to test the theory, we perform event driven 
molecular simulations of a two-dimensional gas of inelastic hard disks, driven by contact with a 
heat bath. Even for strong dissipation, high densities, and small numbers of particles, we find 
that continuum theory describes the system well. With a bath that heats the gas homogeneously, 
strong velocity correlations produce a slightly smaller energy loss due to inelastic collisions than 
that predicted by kinetic theory. With an inhomogeneous heat bath, thermal or velocity gradients 
are induced. Determination of the resulting fluxes allows calculation of the thermal conductivity 
and shear viscosity, which are compared to the predictions of grarmlar kinetic theory, and which can 
be used in continuum modeling of granular flows. The shear viscosity is close to the prediction of 
kinetic theory, while the thermal conductivity can be overestimated by a factor of 2; in each case, 
transport is lowered with increasing inelasticity. 



I. INTRODUCTION 

Collections of inelastically colliding macroscopic par- 
ticles, when driven sufficiently by an external force, ex- 
hibit fluid like behavior such as flow [|l]j2| and instabil- 
ities Although only a small number of particles 
(10^ — 10^) are typically involved in a flow, continuum 
methods are popular tools in modeling, because a cen- 
tury of fluid dynamics experience can be brought to bear 
on the problem. With continuum equations velocity pro- 
files and transfer rates can be calculated, and stability 
analyses performed. Often, plausible continuum equa- 
tions are simply posited, but these equations typically 
contain unmeasurable free parameters. A more rigorous 
approach derives continuum equations from the kinetic 
theory of dissipative gases In principle, a closed 
system of partial differential equations results, analogous 
to the Navier-Stokes equations, with all transport coeffi- 
cients given. 

However, the continuum equations for granular media 
do not share the stature of their molecular fluid analogs. 
For both grains and molecules, the general form of the 
equations is not in question, but the constitutive equa- 
tions relating fluxes of momentum and energy to gradi- 
ents may differ from the simple cases of Newton's vis- 
cosity law and Fourier's heat law. The transport coef- 
ficients of liquids are not routinely calculated from ki- 
netic theory, but are measured experimentally; such mea- 
surements have not been carried out for granular media, 
where the granular temperature (the kinetic energy as- 
sociated with the fluctuational velocities of particles) is 
difficult to control and difficult to measure. In addition, 
the small number of particles and the dissipative nature 
of the medium have led some researchers to the conclu- 
sion that continuum approaches are doomed to failure M . 



As a result, the central question remains open: Can 
continuum equations, derived from kinetic theory and 
supplemented by measurements, model rapid granular 
flows to the same level that Navier-Stokes equations 
model the flows of liquids? Molecular dynamics simu- 
lations will play a crucial role in answering this question, 
since the data they provide can be used to quantitatively 
test the assumptions and the results of kinetic theory. 

The main assumptions, which we will elucidate in 
Sec. H, are: 

1. Single particle distribution functions are nearly 
Boltzmann, 

2. Molecular chaos — Particle velocities are uncorre- 
lated, and 

3. Particle positions are correlated in accord with 
Eq. (^) , the Carnahan and Starling relation for elas- 
tic particles. 

The main results, also discussed further in Sec. ||, are: 

1. The equation of state (Eq. (^3|)), 

2. The constitutive relations (Newton's stress law and 
Fourier's heat law — Eqs. (11) and (|l^), and 



3. The values of the shear viscosity /i, the thermal con- 
ductivity K, and the loss rate of granular tempera- 
ture due to inelastic coUisions, 7 (Eqs. (p^, (p^), 
and (III ' ~ 
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We will use molecular dynamics simulations to test 
these points, and to measure the transport coefficients to 
be used in continuum analyses of granular flows. We have 
found that the simulations quantitatively reproduce ex- 
periments on standing waves in oscillated granular media, 
producing the correct wave patterns and wavelengths |^ 
and secondary instabilities With the numerical sim- 
ulation, we not only have access to experimentally un- 
measurable quantities ||l0[| , but we also can study systems 
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that are not experimentally realizable, allowing us to test 
both the assumptions and results of the granular kinetic 
theory. Once the constitutive relations have been tested 
and possibly enhanced through the use of particle simula- 
tions, direct comparison between continuum theory and 
laboratory experiment becomes possible. 

Laboratory experiments of granular kinetic theory 
have not yet proceeded beyond measuring the single 
particle velocity distribution function [pl]-p^. Simula- 
tions that make contact with kinetic theory have fo- 
cussed mainly on the homogeneous cooling state, a time- 
dependent state that eventually becomes spatially in- 
homogeneous. In those simulations, long-range velocity 
correlations develop and molecular chaos no longer 
holds. A more sophisticated form of kinetic theory, ring 
kinetic theory, in which correlations in particle velo city 
are accounted for, is required for a full description [[l5[ . 
Finally, a number of simulations have been performed on 
one-dimensional granular media in contact with a heat 
bath 1 16-1^] to produce a steady state. In one-dimension 
(ID), spatial IQ and velocity correlations can de- 
velop, and the single particle velocity distributions may 
be nongaussian |18| . The present work extends these ID 
steady state and 2D decaying simulations by examining 
2D steady states with the goal of quantitatively testing 
the assumptions and results of granular kinetic theory. 

In Section III we describe the molecular dynamics sim- 
ulations and the forcing methods. Section [V discusses 



simulations of granular media in which the heat bath is 
spatially homogeneous. In that section, we will check 
the assumptions about the nature of the single particle 
distribution function and the correlations of position and 
velocity, as well as measuring the equation of state and 7. 
In Section ^ we turn to our main purpose: allowing the 
heat bath to vary spatially so that steady state inhomoge- 
neous states may be induced and described. From these 
simulations we can study the constitutive relations, the 
shear viscosity, and the thermal conductivity. Section ^ 
contains concluding remarks. 



II. KINETIC THEORY 

We begin with a brief review of the kinetic theory of 
granular media, which differs only slightly from the ki- 
netic theory of elastic particles as presented in textbooks 
such as p9[ |. The number of particles in a volume and 
velocity element drdv centered at position r and velocity 
V IS given by /(i)(r,v,i)dr(iv; /(i)(r,v,t) is called the 
single particle distribution function. Continuum quanti- 
ties are given as averages over f^^^^ (r, v, t). In particular, 
the number density, average velocity and granular tem- 
perature are defined respectively as 



n{r,t) = J /(i)(r,v,t)dv, 

Mr,t) = - f /(iHr,v,t)vdv, 
n J 



(1) 
(2) 



(3) 



where D is the number of dimensions. Note that the 
granular temperature T is not the thermodynamic tem- 
perature of the particles due to the random motions of 
their molecules, but the analogous kinetic energy due to 
the random motions of the macroscopic particles them- 
selves. 

The evolution of /^^^(r,v,t) depends on the joint 
probability distribution, /^^^(ri, vi, r2, V2, i); collisions 
of two particles change the single particle distribution. 
As in molecular kinetic theory, Boltzmann's assumption 
of molecular chaos is made, i.e., velocities are assumed 
to be uncorrelated, although for sufficiently dense media, 
positional correlations are allowed: 



/^"^(ri,vi,r2,V2,i) = 
<7(a,j.)/(i)(r2-ak,vi,t)/(i)(r2,V2,i), 



(4) 



where a is the particle radius, k is a unit vector pointing 
from the center of particle 1 to the center of particle 2, 
and v = jmra"^ is the solid fraction in 2D. 

The positional correlations are accounted for through 
g(r, v), the radial distribution function, which is defined 
as the probability of having a pair of particles whose rel- 
ative distance lies in the interval r,r + dr, normalized 
by the probability for an ideal gas. This function, evalu- 
ated at the point of contact r = a, gives the increase in 
the probability of collisions due to dense gas (excluded 
volume) corrections. For elastic hard disks, spatial cor- 
relations are described by the formula of Carnahan and 
Starling (20I, 



where 



G{v) = vg{a,v). 



(5) 



(6) 



Equation (^ works well for elastic particles with solid 
fractions below 0.675, where a phase transition takes 
place ]2l|], and is often used in modeling granular me- 
dia 0]. Equation (^) is the definition of G in terms of 
the (unknown) radial distribution function g[r,v), eval- 
uated at r = CT, while Eq. (^ is a particular model for 
G, denoted by the subscript CS. 

An unforced collection of molecules approaches a 
Boltzmann distribution, 



/(i)(i-,v,i) = n(27rr)-^/2g- 



C^/2T 



(7) 



where C = |v— Vo|. Away from equilibrium, the local dis- 
tribution for elastic particles is nearly Boltzmann. Gran- 
ular media dissipate energy with each collision, so that 
the equilibrium state of an unforced granular medium is 
that of no relative motion. For grains, the single particle 
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distribution function is simply assumed to be nearly a 
Boltzmann distribution. 

With these assumptions, and the additional assump- 
tion that the coefRcient of restitution e is only slightly 
less than 1 (particles are only slightly inelastic), equa- 
tions for the continuum mass, momentum, and energy 
can be derived for disks in two dimensions MM: 



— -f V • nvo = 
at 

n—^ + nvo • Vvo -V • P 

ot 

dT 

n— + nvo • VT = -V • q - P : E - 7, 

OT 



(8) 
(9) 
(10) 



where Eij = ^{diVoj + djVd) are the elements of the 
symmetrized velocity gradient tensor E. The constitutive 
relations for the pressure tensor P and heat flux q are 



and 



P = (P - 2ATrE)I - 2/i(E - (TrE)I) 



(11) 



(12) 



where Tr denotes trace and I is the unit tensor. The 2D 
equations close Q with the equation of state, which is 
the ideal gas equation of state with a term that includes 
dense gas and inelastic effects, 



P = {A/Tra^)vT[l + (1 + e)G{v)], 



(13) 



and the predicted values, denoted with a subscript 0, for 
the bulk viscosity, A, 



Ao 



the shear viscosity, /i. 



8iyG{v) It 



TTa 



(14) 



/^o - f + 2 + (1 + -)G(^)] J-, (15) 

2cr G(i^) TT V TT 

the thermal conductivity, k, 

27/1 Q 4 Ft 

^0 = -[7^ + 3 + (- + -)G{^)]J-, (16) 

and the temperature loss rate per unit volume, 7 

3/2 



70 = l-e2 



(17) 



Because of the assumption of near elasticity, the coef- 
ficient of restitution enters only in the equation of state 
and in the expression for the temperature loss due to in- 
elastic collisions. To this order, the thermal conductivity 
and shear viscosity are the same as those given by the 
Enskog procedure for elastic disks 



III. DRIVEN GRANULAR MEDIA 
SIMULATIONS 

We perform event-driven simulations ||2^,Q of a gran- 
ular gas in a two-dimensional periodic box of side length 
L = 52.6(T, in contact with a thermal bath that stochas- 
tically heats particles throughout the volume. Between 
collisions, particles travel freely. In our model the col- 
lisions are instantaneous and binary; they conserve mo- 
mentum and dissipate energy. Particles are assumed to 
be frictionless; particle friction can be incorporated into 
kinetic theories |p5| , p6| , and was included in the simula- 
tions of oscillated granular media, |^,|l^,^ , but introduces 
complications that we wish to avoid for this study. 

When particles collide, new velocities are calculated 
by reversing the component of the relative particle veloc- 
ity along the line joining particle centers and multiplying 
it by the coefficient of restitution e, which is between 
and 1. If e is independent of collision velocity, a fi- 
nite time singularity can occur in the collision frequency, 
a phenomenon known as inelastic collapse |^ , ^ . To 
avoid this simulation-ending occurrence in a convenient 
and natural way, we allow e to vary as 



e(w„) 



1 - Bv^l , Vn < Va 
e ,Vn> Va 



(18) 



where Vn is the component of relative velocity along the 
line joining particle centers (normal to the contact sur- 
face), B — {1 — e){va)~'^ , P = 3/4 and e is a constant, 
chosen to be 0.7. In simulations of oscillated granular 
media, the results are not sensitively dependent on Va 
or /?, and e = 0.7 produces good agreement with ex- 
periment |4p|,p^. In addition to forestalling collapse, 
variation of e allows us to further probe granular kinetic 
theory, which assumes that e « 1. As T varies, the rela- 
tive collision velocity will vary; hence so will e, and the 
relative importance of that assumption can be gauged. 
At a given temperature, a distribution of collision veloc- 
ities and a corresponding distribution of coefficients of 
restitution occur. Varying T varies not only the average 
value of e, but also the amount of deviation around that 
average. 

The variation in e gives rise to a velocity scale that is 
not present in the elastic case. All quantities given below 
are nondimensionalized with the particle diameter a and 
the crossover velocity Va at which the coefficient of resti- 
tution becomes a constant. In particular, the granular 
temperatures all scale with u^. 

Because inelastic collisions remove energy from the sys- 
tem, we must constantly add energy to achieve any sort 
of steady state. The situation is opposite that in sim- 
ulations on nonequilibrium systems of elastic particles, 
where the constant energy input from the driving must 
be removed through an artificial means |2^ . 

Stochastic heating is performed in one of three ways: 
A) White Noise — random kicks (5v are added to par- 
ticles' velocities, B) Random Accelerations — particles 
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accelerate between collisions, and C) Boltzmann Bath — 
particle velocities are obliterated and and replaced with 
velocities chosen from a Gaussian distribution. We dis- 
cuss the motivation and implementation of each in turn. 

A. White Noise 

Williams and Mackintosh |^ introduced white noise 
as a thermal bath for dissipative granular media. In their 
model, a random velocity is added to each particle's ve- 
locity during each time step At. The velocities added to 
each particle are not correlated with one another, nor are 
they correlated with the velocities added in the previous 
At. This model has the advantage that the equation of 
motion for particles between collisions may be written 
down as the Langevin equation 



where Xi is the position of the z-th particle and Q is a 
Gaussian white noise term, i.e., {Ci{t)Q{'t')) — "^FSijSlt— 
t'). The fact that the heating can be analytically ex- 
pressed makes its inclusion into kinetic theory possible. 

This forcing is straightforward to include in our simula- 
tions. Rather than adding the random kicks to particles 
all at once, we kick 2rk randomly chosen particles every 
time there is a collision. The number r^, then, repre- 
sents the ratio between the rate of kicks and the rate of 
collisions. If the kicks are totally random, the center of 
mass momentum of the system will fluctuate, but we de- 
sire that the heat bath can only change the fluctuational 
velocity, not the mean. To ensure that the mean velocity 
remains fixed, we apply random kicks to particles: 

V, ^ V, + |5v|f,,l < i < rfc. (20) 

The kicks themselves are all of the same size, |^v|, but 
the directions f ; are randomly chosen. Then, kicks in the 
opposite directions are applied to another rk randomly 
chosen particles: 

Vt v., - |5v|fi_.rfc,r-fc < i < 2rfe. (21) 

Because each kick requires recalculation of the kicked 
particle's collision list, we want to minimize the kick fre- 
quency. Empirically, we find that our results are inde- 
pendent of rfe for Tfc > 1 and rfc((5v)^ constant. This is 
not surprising, since the average length of randomly 
oriented kicks of length 6v is y^(5u, and only the to- 
tal kick between collisions matters when collisions occur. 
For this reason, we perform most of our simulations with 
Tfc = 1, where the collision rate equals the kick rate. 

B. Random Accelerations: The Air Table Model 

While white noise forcing has the advantage that it can 
be incorporated into the kinetic theory relatively easily, 



it has the disadvantage that it does not model any partic- 
ular real system. In most experiments, energy is added 
to the granular media through a boundary, causing gra- 
dients in the energy perpendicular to that boundary. For 
a vertically oscillating layer in a gravitational field 
the situation is even more severe; although one might 
like to imagine the oscillating wall as thermal, provid- 
ing stochastic kicks to the layer, the interactions of the 
layer with the plate are strongly correlated to the plate's 
motion. The plate introduces both spatial and temporal 
gradients. 

One experimental system is capable of producing ho- 
mogeneous steady states, specifically, a collection of 
pucks on an air table [ p^||3^ . Although the purpose of 
the air is to levitate the pucks, it also drives their hori- 
zontal motion. Either due to inhomogeneities in the air 
fiow or because the pucks are not perfectly parallel to the 
surface of the table, the pucks accelerate uniformly from 
one collision to another ||l3|. This driving counters the 
loss of energy due to inelastic collisions and produces a 
steady state. 

We model this air table driving by allowing each par- 
ticle to move under a uniform acceleration: 

a-i = aor, (22) 

The magnitudes of all particle accelerations, oq, are the 
same, but the directions, f^, are randomly and uniformly 
chosen. The direction of a particle's acceleration changes 
stochastically. When a collision occurs, 2rk particles are 
given new f^. In order to conserve total momentum, we 
hold the total acceleration of the particles at zero by giv- 
ing exactly opposite accelerations to pairs of particles. 
Initially, each particle is paired with another, and these 
are given opposite accelerations. Later, when one parti- 
cle is chosen and its acceleration randomized, its partner 
particle is also given a new acceleration, opposite to the 
first particle's. 

Experiments suggest that particles accelerate uni- 
formly from collision to collision, with most of the 
changes in acceleration happening at collisions 
Therefore, = 1 is probably a relatively good model 
for the air table experiments. As rk increases, the par- 
ticles feel a constant acceleration over a small temporal 
range. In the limit that — > oo, the model is the same 
as the white noise model, which is rate independent. In 
practice, we observe that the distribution of collision ve- 
locities for the accelerated particle model approaches that 
for the white noise model at « 8. 

Variation of coefficient of restitution with normal col- 
lision velocity is critical for simulations with the accel- 
erated forcing, since inelastic collapse-like collision se- 
quences are more prevalent. After a collision, particles 
move away from one another. Without relative acceler- 
ation, they must collide with other particles to reverse 
their velocity if they are to recollide. Particles with rela- 
tive acceleration, however, may collide and recollide with- 
out striking another particle in the intervening time, just 
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as a ball bounces repeatedly on the ground. On average, 
the relative acceleration is likely to change during this 
interval, but that is not guaranteed for any given pair. 
Further, since accelerations are changed when collisions 
occur, the rate of acceleration fluctuation is dependent 
on the collision rate in the entire medium. A given pair 
of particles can hijack the collision sequence of the gas, 
rapidly recoUiding with one another. If the coefficients 
of restitution are constant, this scenario will produce col- 
lapse. By allowing the collisions to become more elastic 
for decreasing relative velocities, however, collapse is pre- 
vented; eventually, the relative acceleration will change, 
and the particles will move apart. 



C. The Boltzmann Bath 

Finally, we introduce a heat bath that approximates 
the assumption of molecular chaos. Molecular chaos as- 
sumes that the velocities of colliding particles are uncor- 
related; particles collide, and before they collide with one 
another again, they collide with a large number of other 
particles, losing the memory of the initial condition. A 
strong heat bath can perform the same function if it re- 
places the particle velocities with new velocities chosen 
from a given distribution. If the heat bath interacts with 
the particles often enough, molecular chaos will be guar- 
anteed. While this bath is wholly unphysical, it produces 
a situation in which the kinetic theory is expected to ap- 
ply exactly; it is a useful check on calculations, and helps 
to elucidate the role of velocity correlations. 

Implementation of the Boltzmann bath is simple. 
When two particles collide, 2rk particles are randomly 
chosen and given velocities chosen from a Boltzmann dis- 
tribution with a specified temperature. 
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FIG. 1. Average coefficient of restitution vs T for f — 0.5, 
for the three forcing methods. +: White noise forcing, <>: 
Accelerated forcing, *: Boltzmann Bath. As discussed in the 
text, T is nondimensionalized with v^- 

The average coefficients of restitution for a number of 
runs at different temperatures are shown in Fig. ^. Note 
that at high temperatures, the average coefficient of resti- 
tution begins to rise for the accelerated forcing. This is 
due to the inelastic collapse-like collision sequence de- 
scribed earlier for this forcing. As the magnitude of 
acceleration increases to produce higher temperatures, 
these multiple collisions become more and more impor- 
tant, leading to a large number of collisions with very low 
velocities, and so a high average coefficient of restitution. 

In the state produced by homogeneous forcing, we 
can measure the single particle distribution function, the 
temperature produced by the forcing, the pressure, the 
radial distribution function, velocity correlations, and 
loss rates. These quantities can be compared to the cor- 
responding quantities for elastic simulations and to their 
assumed or calculated values from kinetic theory. 



IV. HOMOGENEOUS FORCING, 
CORRELATIONS, AND THE EQUATION OF 
STATE 

With any of the thermal baths just described, we can 
set up inhomogeneous states by varying the strength or 
rate of forcing over space. However, in the simplest case 
we force homogeneously. The simulations are performed 
for a variety of solid fractions, forcing rates, and forcing 
strengths for the three thermal baths. 



A. Single Particle Velocity Distributions 

The lowest order approximation to / in the kinetic 
theories is usually chosen to be a Boltzmann distribu- 
tion, the form of / for an undriven elastic gas. A driven 
inelastic gas, although it approaches a steady state, is 
by no means guaranteed to act like an undriven elastic 
gas. Nevertheless, the single particle distribution func- 
tions measured from the simulation are all close to Boltz- 
mann distributions, as seen in Fig. H 
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T = 1 .05, V varies v = 0.5, T varies 
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FIG. 2. Single particle distribution function Q{v') for (a) 
and (b) White noise forcing, (c) and (d) Accelerated forcing, 
and (e) and (f) the Boltzmann bath, all with = 1. The ve- 
locities are scaled with the temperature T, so that v' = v/VT 
and Q{v') — Pr(v')\/W, where Pr{v') is the probability dis- 
tribution of v' . In the left column, the average temperature 
is approximately 1.05, and the solid fraction is varied 
u = 0.1, *: u = 0.4, o: V = 0.6, A: u = 0.8.) In the 
right column, i/ is fixed at 0.5 and the temperature is var- 
ied; (b) r = (-f)1.93 X 10~^ (*)3.13 X 10-^ (o)1.06, (A)1067. 
(d) r = {+)3.0 X 10"^(*)l.l X 10"^(o)1.05, (A)256. (f) 
T = {+)1.2 X 10~^ (*)1.1 X 10"^ (o)1.02, (A)102. The solid 
curves are Boltzmann distributions. 
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v' 

FIG. 3. The velocity distribution function Pr of v' = v/VT 
from a simulation with accelerated forcing at u — 0.5 and 
T — 1.05, divided by PruB, a Maxwell-Boltzmann distribu- 
tion with T = 1.05. The two curves are for the two velocity 
components. 

Overall, the accelerated forcing produces the strongest 
deviations from Maxwellian, and the Boltzmann bath, 
unsurprisingly, produces the least deviation. In all cases. 



the deviations become stronger as the density and tem- 
perature increase (recall that increasing temperature has 
the same effect as decreasing the average coefficient of 
restitution). These deviations tend to flatten the distri- 
bution, increasing the probability in the tails and slightly 
in the peak, and decreasing the probability in between, 
as displayed in Fig ^. Similar types of deviations, but 
much stronger, have been observed in experiments on a 
dilute, vertically oscillated granular layer [ pT[ . 

B. Equation of State and the Radial Distribution 
Function 

The equation of state, (p^, relates the pressure, den- 
sity and temperature to the coefficient of restitution and 
G{i'). The virial theorem of mechanics as applied to hard 
spheres can be used to calculate the equation of state 

Py = 7VT+-^^k.Av„ (23) 

where the sum is over all collisions that occur during the 
measurement time tm, Av^ is the change in the velocity 
of the i-th particle due to the collision, and k is the unit 
vector pointing from particle center to particle center. 
In this form, measurement of pressure reduces to mea- 
surement of the average particle energy and the average 
change in the normal velocity at collision; we measure 
pressure with this method. 

Using Eq. ( p3| ) to measure pressure, and assuming the 
equation of state (Eq. (p^), we produce a measurement 
of G'(z^), denoted Gs(i^), where the subscript s stands for 
simulation. This measured value of G will be compared 
to the Carnahan and Starling value Gcs{^) from Eq. (j|). 
Accurate characterization of G{i') is important, because 
it occurs in the expressions for transport coefficients. 
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FIG. 4. (a) Ga{v) for inelastic hard discs driven by +: 
white noise forcing, o: accelerations, and *: Boltzmann forc- 
ing. The solid curve is the Carnahan and Starling relation 
Gcs{v), given by (|). (b) The ratio of G,{v) to Gcs{v)- AH 
runs have = 1 and T = 1.05. 

We calculate Gs{v) for the three types of forcing as v 
varies. The results are shown in Fig. IJ. For v below w 
0.675, where elastic particles undergo a phase transition 
to an ordered state |2^, the white noise and accelerated 
runs produce lower G than elastic runs; Boltzmann runs 
have Gs{v) slightly above the elastic values. 

As the temperature decreases, e 1, and the values 
of Gs(j^) must approach the elastic values. Therefore, 
Gsiiy) must be temperature dependent; this dependence 
is shown in Fig |^, along with the value of Gs(i^) given by 
Eq. (||). As T decreases, the inelastic Gsii') approaches 
the elastic G, and at high T, where e is independent of 
T, Gs(y) becomes independent of T. As for the single 
particle distributions, the accelerated forcing shows the 
greatest deviation from the elastic behavior. 
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FIG. 5. Gsiy) vs T for u = 0.5. +: white noise, o: accel- 
erated, *: Boltzmann. The dotted line is the Carnahan and 
Starling relation Gcs{v), given by (jH]), for v — 0.5. 

The substantial differences between Gs(z^), deduced 
from the equation of state, and Gcsii^Yt predicted by 
the Carnahan and Starling relation (Eq. (^)), comes from 
two sources. First, Eq. (|^) may not give the correct value 
of the radial distribution function at the point of con- 
tact, f/(cr, v), for granular media. Second, Gs{v) may not 
equal yg{(J, v). Because Eq. is the definition of G, this 
would amount to a change in the equation of state, which 
we assumed to be correct in the calculation of Gs{v)- 
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FIG. 6. g{r) a.t u = 0.5, T = 1.05 for (a) white noise 
(rfc = 1) and (b) Boltzmann (r^ — 4) forcing. The dot-dashed 
lines represent the value of g given by the Carnahan and Star- 
ling relation Eq. ^ for v = 0.5, while the dashed lines show 
Gsi^)/!/. For white noise forcing, g{(T,v) coincides with nei- 
ther line, while for the Boltzmann bath, g{a, v) coincides with 
Gs{v)/u. 

To elucidate these two causes, we plot g(r, v) for a 
run with white noise forcing and a run forced with a 
Boltzmann bath. Each plot also indicates the value of 
g{i7,v) predicted by Eq. (H), as well as that predicted 
l3y Eq. (^, assuming G{v) — Gs{v)- For neither forc- 
ing type does (^), the Carnahan and Starling relation for 
Gcs{v), properly predict g{a,v)] rather, inelastic parti- 
cles are more likely to be nearly in contact than elastic 
particles at the same density and temperature. Further- 
more, while Gs{v) = vg{(T, v) for the Boltzmann forcing, 
this does not hold for the white noise forcing or the accel- 
erated forcing (not shown). Even though g{a^ v) is larger 
than that predicted by the Carnahan and Starling rela- 
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tion, G's(i/) < Gcs{>^), indicating that the equation of 
state is incomplete. Recalhng that the Bohzmann driv- 
ing represents particles in contact with a highly random- 
izing bath, we conclude that the failure of the equation of 
state, dl^), is due to incomplete randomization of parti- 
cle velocities through collisions, or in short, a breakdown 
of molecular chaos. 



C. Velocity Correlations 

Molecular chaos is the assumption that particle veloc- 
ities are uncorrelated. Knowing the velocity of one of a 
pair of a colliding particle gives no information about 
the velocity of the other. In light of the behavior of 
G, and simulations of driven ID and cooling 2D gases 
that showed strong velocity correlations , we measure 
velocity-velocity correlation functions. 

Given two particles, labeled 1 and 2, k is a unit vector 
pointing from the center of 1 to the center of 2. Particle 
I's velocity then has a components w[' parallel to and Vi 
perpendicular to k; likewise for particle 2. We define two 
correlation functions 



Y^vivi/Nr 



(24) 
(25) 



where the sums are over iV^ particles such that the dis- 
tance between the two particles is within 5r of r. If par- 
ticle velocities are uncorrelated, ^^^'^ '^^^1 
both give zero. 



FIG. 7. Velocity correlations as a function of particle sepa- 
ration at u — 0.5, T — 1.05 for: white noise, (o) accelera- 
tions, (x) Boltzmann forcing, and (A) elastic particles. Each 
curve is built from around 100 frames separated in time by 
100 collisions per particle, and 5r = cr/10. Both the elastic 
particles and the particles forced with the Boltzmann bath 
have essentially zero correlation over most of the range. The 
Boltzmann bath shows positive correlations only at very short 
range. 

The parallel and perpendicular velocity correlations 
are plotted in Fig. |7 for the three types of forcing and 
for elastic particles. Both for particles driven with white 
noise and accelerations, strong long-range velocity corre- 
lations are apparent, with more correlations produced by 
the accelerated forcing, consistent with its stronger devi- 
ations in the single particle velocity distribution and in 
G. These correlations are not small, reaching as much as 
40% of the temperature; typically, the perpendicular cor- 
relations are about one-half of the parallel correlations. 
Further, these correlations are long range — they extend 
the full length of the system. The parallel correlations 
drop to zero at L/2, while the perpendicular correlations 
reach zero around r ~ lOcr, and have a negative value 
but zero derivative at L/2. The long-range of nature of 
the correlation is not due to the size of the computational 
cell. Similar cell-filling correlations were observed in runs 
4, 16, and 64 times larger [3^ ]. 

For the Boltzmann forcing, some correlations are visi- 
ble at very short range; inelastic collisions are trying to 
establish correlations, but before these correlations can 
be extended to long range they are wiped out by the 
thermal bath. 
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D. Loss Rate 

The loss rate of temperature due to inelastic collisions, 
7, divided by the rate calculated from kinetic theory, 70 
(see Eq. (17)), is shown for the three forcing methods 
as a functions of T in Fig. ||(a) and i' in Fig. ||(a). For 
the calculation of 70, G was taken from the equation 
of state measurements, and the average e within a run 
was used. Most surprising is the increased loss rate over 
kinetic theory at relatively low temperature. 
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FIG. 8. Temperature dependence of loss rate with v = 0.5. 
(a) 7/70, where the kinetic theory result 70 [Eq. (p7|)] assumes 
a velocity-independent restitution coefficient, (b) 7/7e, where 
7e takes into account the velocity dependence of e. The sym- 
bols denote the type of forcing: -I- White noise; <> Accelerated; 
X Boltzmann. The dotted lines show ■j^{v^)c/{vn)c for the 
white noise and accelerated runs; see (p2|). 

The calculation of 7 requires only the evaluation of 

a r r 1 
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-{1 - e2)(vi2 • k)3/^'Hvi, V2)c;Mvidv2, 



(26) 



where dfl is the angle element. This expression simply 
averages the energy lost per collision, 



e>^ = i(l-e2)K2-kf (27) 



over all possible collisions. The remaining factor of 
(vi2 • k) takes into account the fact that grains travel- 
ing towards one another more rapidly are more likely to 
collide during any given time interval. 

The kinetic theory result for the loss rate, Eq. (17), 
follows from the coUisional integral, Eq. (p6|), under two 
assumptions: molecular chaos, which ought to be a more 
reasonable assumption at lower temperatures, and the 
independence of e on the variables of integration, so that 
it is pulled out of the integral as a constant. At lower 
temperature, where the variation of e with Vn leads to 
a distribution of e at a given temperature, the second 
assumption fails. As temperature drops still farther, the 
spread of e reduces, since e is bounded from above by 1; 
at the very lowest T, 7 does approach 79. 



Substituting the function form of e(t;„), Eq. (|T^), into 
the collisional integral, Eq. (|2^), assuming molecular 
chaos, and performing the integrations, we arrive at an 
equation for 7 = 7e that takes into account the variation 
of e with Vr, : 



^■3^3/2 



[(1 - + AT) eM-vl/'iT) + 4/] , 

(28) 



where 



- A'^2^f^T^+'^{r{2 + /3) - r(2 + P, vl/4T)), (29) 

r(a) is the gamma function, and r(a, b) is the incomplete 
gamma function. In the limit that Va ^ 0, 7e — > 70- 
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FIG. 9. (a) 7/70 and (b) 7/7^ versus u a.t T = 1.05 The 
symbols denote the type of forcing: + White noise; o Accel- 
erated; X Boltzmann 

Figs, ^(b) and Wh) show 7/7e for the same simula- 
tions shown in FigsT^(a) and |9|(a). Taking the variations 
in e into account removes the underestimation of 7. In 
the revised picture, 7 approaches 7e at low T and at low 
v, but as either increases, 7 drops from the value pre- 
dicted by 7e . This is due to the velocity correlations pro- 
duced by the inelasticity. Locally, particles are moving 
together, reducing collision velocities and collision fre- 
quencies, thereby reducing the loss rate; see Fig |l^. For 
the Boltzmann forcing, velocity correlations are wiped 
out, and 7 is close to 7e. 
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V. INHOMOGENEOUS FORCING AND 
TRANSPORT COEFFICIENTS 

So far, we have been concerned with homogeneous forc- 
ing. By using appHed forcing that varies spatially, we 
can induce inhomogeneous steady states. Then, by mea- 
suring fluxes, transport coefficients are calculated, and 
compared to kinetic theory. Inhomogeneous states have 
only been calculated for the accelerated forcings; mea- 
surements for the homogeneous state show that devia- 
tions from the elastic case are strongest in this case. 



FIG. 10. Probability distribution of Vc = |vi — V2I, the 
magnitude of relative velocities at collision for different forc- 
ings white noise, (o) acceleration, (x) Boltzmann, and 
(A) for elastic particles, all at u — 0.5 and T — 1.05. The 
solid curve is the distribution predicted by uncorrelated colli- 
sions between particles chosen from Boltzmann distributions. 
Positive velocity correlation of nearby 



_ /4T 



particles causes a reduction in the collision velocities, 
hence a reduction in 7. 



and 



Writing 7 in terms of average quantities makes its de- 
pendence on the collision velocity more explicit. The loss 
rate is identically equal to the average energy lost per col- 
lision, {AE)c, times the average collision frequency per 
volume f/V 



1 



^ = {AE)J/V^-{l-e'){vt,)J/V, 



(30) 



where the c subscript denotes an average over collisions, 
and assuming again that e is independent of collision ve- 
locity. Similarly, the virial equation of state, Eq. (|23|), in 
terms of average quantities is 



P = {4/Tra'^)vT{l ■ 



V(J 1 + 
d ■' 2 



(31) 



Solving for / in terms of G from Eqs. ( p^ ) and (|3^), and 
substituting this into Eq. (^ ) we obtain 
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nT. 



(32) 



If the distribution of relative normal velocity at col- 
lision is equal to that predicted by molecular chaos, 
P(w„) = (l/2T)^;„exp-w2/4T, then {vl)c = 4T and 
{vn)c = VttT, so that 7 = 7o is recovered. 

As seen in Fig. ^ however, the distribution of collision 
velocities is different from the molecular chaos values due 
to the presence of velocity correlations. To show that this 
deviation accounts for the remaining difference between 
7 and 7e, (i'^)c and {vn)c were calculated in the simu- 
lations. Their ratio, normalized by the molecular chaos 
value TT is plotted on Figs. ||(b) and ^(b) as dotted 

lines. Except where there are wide variations of e within 
a given run (such as at high velocity in the accelerated 
runs) , the change in the relative collision velocities tracks 
the change in 7. 



A. Thermal Conductivity 

Recall that with the accelerated forcing, the direction 
of the accelerations of particles fluctuated at a fixed rate, 
but the magnitude was always the same. To induce a 
thermal gradient in the simulation, we allow the mag- 
nitude of the acceleration to vary in space. Specifically, 
when the acceleration of a particle is to be rotated, the 
magnitude of its acceleration is given by 



«o(l-|(^-^i)/(^i)|), 



(33) 



i.e., we apply a linear gradient in the forcing, dependent 
on one spatial direction (z), peaked in the center of the 
cell and falling to zero at the periodic boundary. In order 
to preserve the center of mass momentum, the partner 
particle receives the opposite acceleration, reg ardle ss of 
its position in the cell, as described in Section IIIB. 



Under this forcing, a stable thermal profile develops, 
as seen in Fig. The system reaches a mechanical equi- 
librium; the pressure is nearly constant in z. Constant P 
and varying T imply varying v through the equation of 
state, and this variation is observed. 
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FIG. 11. Profiles of *: T, □: u, and <>: P/5 for inhomo- 
geneous forcing in which the magnitude of the acceleration 
depends linearly on the distance from z/L — ^. The pressure 
is nearly constant, so that the variation in T induces a vari- 
ation in u through the equation of state. For this run, the 
average solid fraction is 0.75. 
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The cell is divided into 50 slabs along z for measure- 
ment purposes; for each slab, snapshots allow calculation 
of the average T and v within. Because the pressure 
may in principle vary over space, its measurement using 
the virial fails. Pressure is measured at the interfaces 
between the slabs by keeping track of the normal (z) 
momentum flux through these boundaries, both due to 
particles freely traveling through them, and due to col- 
lisions between particles that lie in different slabs. In 
addition, the energy added due to the forcing and the 
energy lost due to inelastic collisions are separately ac- 
counted for each slab. The difference between the energy 
gain and loss in a given slab must, in a steady state, be 
made up for by the difference in energy flux through its 
two boundaries, so that the net rate of change of the en- 
ergy in the slab is zero. Assuming that the energy flux 
through the line at z = (0, L) is zero due to symmetry 
allows calculation of q^, the heat flux through each slab 
boundary; see Fig 
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FIG. 12. The x : thermal gradient and A: heat flux for 
the run shown in Fig. 0. 

Once Qz and dT/dz are calculated, Fourier's heat law, 
Eq. (p^, can be used to calculate the thermal conduc- 
tivity K. The results of many simulations, holding the 
average v at 0.75, but varying oq in Eq. ( p3| ) and there- 
fore the size of the thermal gradient, are shown in Fig. |l^, 
are compared to the result of Enskog theory, as given by 
Eq. dl^). At low temperatures, Enskog theory does a 
good job predicting k. However, as the temperature in- 
creases and e{vn) decreases, Enskog theory does worse; at 
the highest temperatures, Enskog theory overestimates 
the thermal conductivity by a factor of 2. 
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FIG. 13. Ratio of k measured from simulations to kq from 
Enskog theory (Eq. ([l6|)). Each symbol denotes a difi'erent 
run, but for each run, the average solid fraction is 0.75. 

Note that this calculation does not test Fourier's Law; 
rather we assume that Fourier's law is correct, and use it 
to calculate k. Analysis based on closures of the Boltz- 
mann equation predict a term in the heat flux propor- 
tional to the density gradient |3^. If such a term had 
a sizeable magnitude and were ignored, it would cause a 
reduction in the observed n. 



B. Shear Viscosity 

Spatial inhomogeneity in the magnitude of the forc- 
ing led to a stationary inhomogeneous temperature field, 
allowing measurement of heat flux and thermal conduc- 
tivity; spatial inhomogeneity in the mean of the forcing 
leads to a stationary inhomogeneous velocity field, allow- 
ing measurement of the momentum flux and the shear 
viscosity. In particular, particle accelerations are chosen 
according to 



ao(0.01sin {2ttz/L) 



(34) 
(35) 



where i/jq and ipi are numbers chosen randomly from a 
Gaussian distribution with zero mean and standard de- 
viation of 1. 
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FIG. 14. Average velocity in the y direction as a function 
of 2 at T = 0.21, v = 0.6. The sohd line is sinusoidal. 

This forcing produces steady states with velocity, tem- 
perature, density, and stress fields like that shown in 
Fig |l^. The velocity profile is nearly sinusoidal, and 
the temperature, pressure, and solid fraction are essen- 
tially independent of z. In the simulations discussed so 
far, we have only considered the scalar quantity T = 
{(v— {v))'^) /D, where D is the number of dimensions and 
the denote averages over particles. This temperature 
is more generally the trace of the temperature tensor: 



Tij = {{vi - {vi)){vj - (vj))), 



(36) 



where i,j range over the directions, and Vi denotes the 
i-th component of the velocity. In principle, Tyy need 
not equal T^^ if the rate at which fluctuational energy 
is traded between the directions is slower than the rate 
at which it is added anisotropically; such is the case in 
a vertically oscillated granular layer, where vertical fluc- 
tuational energy can be twice that of the horizontal. In- 
troducing a bias in the acceleration of only 1%, however, 
does not introduce anisotropy into T; for larger shear 
rates this is not the case. To study the simplest case, we 
restrict our simulations to biases of 1%. 



In section V A we described calculation of the pressure 



by measuring the normal momentum flux through planes. 
Measuring the tangential flux through the planes, and in- 
troducing a second set of planes orthogonal to the first, 
allows calculation of the full 4-component pressure ten- 
sor. As seen in Fig. |l^, the pressure tensor is anisotropic 
even though T is isotropic; the anisotropy increases as T 
increases, or e decreases. For e w 1, the stress difference 
is approximately proportional to 1 — e, but the variation 
in e within a given run probably plays a role, as it did in 
the loss rate; see Fig. [l5[ 
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FIG. 15. Normal Stress difference divided by pressure as a 
function of 1 — e. The line has slope 1. 

For each run at fixed T, we can test Newton's viscosity 
law, 



P --n^ 



(37) 



where the viscosity /i is a constant of proportionality. 
A typical result is shown in Fig. |l^, where the linear 
relation of Eq. ( ^ ) is shown to hold. The slope of these 
curves then provide values for /i which can be compared 
to the Enskog result from Eq. (15); the results are shown 
in Fig. O. 
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FIG. 16. The yz component of the pressure tensor versus 
the corresponding velocity derivative for the run of Fig. |l^. 
The slope of the best fit line is — /x. 
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FIG. 17. Viscosity, normalized by the Enskog value ^q, as 
a function of T, for v = 0.6. 

Unlike the loss rate and thermal conductivity, we find 
that the Enskog theory underestimates the shear viscos- 
ity at lower temperatures. However, the trend of de- 
creasing transport with increasing T is the same. Even 
for elastic particles, Enskog theory is not expected to 
work to arbitrarily high solid fractions; as density in- 
creases, deviations from Enskog theory are expected. As 
the inelasticity of particles increases, velocity correlations 
increase, reducing coUisional momentum transport and 
lowering the viscosity. 
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VI. CONCLUSION 

Volumetric driving of granular media leads to nearly 
stationary states that are amenable to comparison with 
kinetic theory, allowing us to test the six points of kinetic 
theory listed in Sec. |. Given that many of our simula- 
tions involve coefficients of restitution that are not near 1, 
the general level of agreement with kinetic theory is sur- 
prisingly good, suggesting that continuum approaches to 
dissipative granular media are capable not only of quali- 
tative, but also quantitative descriptions of real systems. 
We now discuss each of the six points in turn. 

Single particle distribution functions are nearly Boltz- 
mann. For all forcing types, temperatures (coefficients of 
restitution), and densities, the single particle distribution 
functions are close to Boltzmann distributions (Fig. |^). 
The deviations from Gaussian (Fig. |3|) are consistent with 
but much smaller than those seen in experiments on thin 
(< 1 layer deep) oscillated granular media [0. In those 
experiments, deviations appear to be due to spontaneous 
spatial variations in temperature that are not taken into 
account in the analysis |35|; such variations become large 
for cooling, unforced granular media. The smaller devia- 
tions exhibited in our simulations may represent smaller 
temperature fluctuations. 

Particle velocities are correlated. Standard kinetic 
theory assumes molecular chaos: particle velocities are 
uncorrelated. In our simulations, as in simulations of 
cooling granular media , strong velocity correlations 
(Fig. ^ with a characteristic vortex structure |33) de- 
velop. In the steady state, these correlations extend over 
the entire cell. Molecular chaos is not required for kinetic 
theory; a closure for the kinetic equation is critical. Simi- 
lar considerations have lead van Noije and Ernst to apply 
ring kinetic theory ||3^ ], which allows for correlations in 
particle velocities, to a granular system 

Spatial correlations are stronger than predicted by the 
Carnahan and Starling relation (Eq. Even for sim- 

ulations using the Boltzmann bath, in which velocity 
correlations are removed, the effect of inelasticity is to 
increase the amount of spatial correlation at r = cr 
(Fig. ^(b)). In other words, particles are more likely in 
the inelastic case to be close to one another than elastic 
particles in the same thermodynamic state. The size of 
this effect is dependent on the inelasticity, but can be 
greater than 15%. 

The equation of state overestimates the collisional con- 
tribution to pressure because it ignores velocity correla- 
tions. The factor in the equation of state describing the 
contribution from collisions, Ggiy), as calculated from 
the measurement of pressure and the equation of state, 
is smaller for white noise and acceleration forcings than 
that predicted by the Carnahan and Starling relation 
(Eqs. (H) and (|)), denoted Gcs{v) (Fig. |). In turn, 
Gcs{v) is smaller than the actual G{v) = vgiv), as dis- 
cussed in the previous paragraph. Because velocity cor- 
relations were ignored in the derivation of the equation 



of state (^_3|), the pressure due to collisions that G de- 
scribes is overestimated. To some degree, the increased 
positional correlation and increased velocity correlation 
work against one another; the first increases the collision 
frequency, while the latter decreases it. The net result is 
that the Gs{v) from the pressure measurement is closer 
to Carnahan and Starling (^ Gnsi^) than if there were 
only velocity correlations (Fig. ^(a)). 

Newton's stress law works well for low stress. Even 
at the highest inelasticity, e = 0.7, no deviations from a 
linear relation between stress and strain rate (Eq. (pl|)) 
were observed (Fig. |l^). However, in order to keep the 
temperature isotropic, we have limited ourselves to cases 
in which < T; many flows of interest are supersonic, 
with average velocities much larger than Vt. Because k 
depends on T, and therefore on position, Fourier's heat 
law was not tested in the same manner that Newton's 
viscosity law was. 

Temperature loss rate, 7, and thermal conductivity, k, 
are reduced by inelasticity, while shear viscosity, /i, is 
predicted relatively well by Enskog theory. For increasing 
inelasticity or density in homogeneously forced runs, ve- 
locity correlations also increase. As a consequence, the 
relative collision velocity decreases (Fig. |l^), leading to a 
reduction in the temperature loss rate due to inelastic col- 
lisions (Figs. 1(b) and |(b)). In the most severe cases ex- 
amined, once corrected for variations in e, this deviation 
could be as high as 20%. Inhomogeneously forced runs 
allowed calculation of /i, and, assuming Fourier's law, k; 
/i never deviated from the prediction of Enskog theory by 
more than 15% (Fig. |l^), while K was found to be smaller 
than predicted by a factor of two for high inelasticities 
(Fig. |l^). This differential success suggests that velocity 
correlations, which should be present in both cases, are 
not responsible for the large reduction in k. Rather, the 
inelasticity itself may be the cause. Enskog theory, ap- 
plied to granular media, assumes that e « 1; the values 
of kq and /io are the same as those for elastic particles. 
When grains collide, energy is dissipated, so that the 
amount transported coUisionally is necessarily reduced 
as e decreases. On the other hand, momentum is still 
conserved, so that ^ is relatively unaffected. Also, some 
deviation from Enskog theory is possible due to a term 
in the heat flux proportional to density gradients. 

The results described above suggest a number of av- 
enues for future research. First, measurements of vis- 
cosity should be extended into the supersonic regime. 
Second, more extensive calculations of thermal conduc- 
tivity at different densities, and with different spatial 
forcings, should be undertaken to ascertain the role of 
density gradients. Third, time-dependent calculations 
should be performed, allowing measurement of the bulk 
viscosity. Such calculations should also provide measure- 
ments of the frequency dependence of the transport co- 
efficients, which may be relevant for oscillated granular 
media. Fourth, the granular continuum equations can be 
used to perform stability calculations on problems such 
as vertically oscillated granular media [p8^,p9|. Finally, 
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new forcing geometries should be explored, allowing di- 
rect comparison between particle simulations, continuum 
theories, and experiments. 
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